home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATHEMAT / 1025.ZIP / LAPLACE.C < prev    next >
C/C++ Source or Header  |  1987-08-15  |  37KB  |  1,207 lines
  1. /*
  2. #define DEBUG1
  3. #define DEBUG2
  4. #define DEBUG4
  5.  
  6.     laplace - solve Laplace's equation for specified boundary conditions
  7.                 in two dimensions
  8.     bugs...
  9.         in mark, a line ending near a grid LINE, not just near a grid point,
  10.         should be registered.
  11.         
  12.     history...
  13.         2 Aug 87    line segments ending near a grid point result in type
  14.                     zero grid points.  type matrix browser installed
  15.                     (Z-100 only).  Several up/down inconsistencies
  16.                     removed.  Crossing display function now shows
  17.                     coordinates of crosspoint.  A type zero grid point
  18.                     is never redefined, even if other crossings are
  19.                     recorded for that point.  Equations for corner
  20.                     points have nonzero coefficients only for neighbors
  21.                     on the grid.
  22.         1 Aug 87    -q switch defeats comments during calculation.
  23.         28 Jun 86    -g switch sets edge of grid to zero potential, 
  24.                     -m switch specifies a margin around the specified boundaries.
  25.         21 Jun 86    Allowing boundary to be between grid points.
  26.         31 May 86    showing boundries and calculated potentials are optional.
  27.  
  28.     approximate performance...
  29.  
  30. h=w=32, n=64
  31. Current time is 20:58:25.76
  32. Current time is 21:00:03.36 -> 5:30.86 or 330.86 sec
  33.  
  34. h=w=64, n=128
  35. Current time is 21:01:11.29
  36. Current time is 21:07:08.87 -> 5:57.58 or 357.58 sec
  37. ...producing 42K file
  38.  
  39.         ...on Z-100 with V20 at 7.5 MHz, all files on ramdisk.
  40.  
  41. h=w=32, n=64
  42. 56 seconds
  43.         ...on IBM AT with 80286 at 8 MHz and 80287 at 5.3 MHz, 
  44.         all files on ramdisk.
  45.  
  46.     notes...
  47.         need to update above performance.
  48.         should investigate recursive subdivision by factors other than 2.
  49.         can't handle boundary crossing edge of grid: referencing nonexistant
  50.             grid points.
  51.         needs option for cylindrical symmetry
  52.         Should optionally output the grid used, in GRAPH input format.
  53.  
  54. */
  55. #define SWITCH
  56.  
  57. #define DESMET            /* deleting this only disables the timing function */
  58.  
  59. #include <stdio.h>
  60. #include <math.h>
  61.  
  62. #define VERSION "1.0"
  63. #define MAX_CROSS 375        /* max # times boundary can cross grid */
  64. #define ENTRIES 375            /* maximum # points in boundary point curves */
  65. #define MAXBOUNDARIES 50    /* maximum # separate curves in boundary */
  66. #define PERMITTIVITY 8.85418e-12    /* permittivity of free space, farad/m */
  67. #define SCALE 16384
  68.  
  69. #define MAX(a, b) ((a)>(b)?(a):(b))
  70. #define streq(s, t) (strcmp(s, t)==0)
  71.  
  72. typedef struct
  73.     {int grid;    /*    cell number */
  74.     int direct;        /*    direction: 0, 1, 2, 3  for  +x, -x, +y, -y  */
  75.     double dist;    /*    distance from cell center to boundary crossing point 
  76.                         (fraction of grid spacing, ranging from 0 to 1)  */
  77.     int volt;        /*    negative of boundary # (so it's an index into volt[]) */
  78.     int dummy[4];    /*    ...so sizeof(crossing) >= sizeof(boundary_rule)    */
  79.                         /* >>> delete? <<< */
  80.     } crossing;
  81.  
  82. typedef struct
  83.     {int var[4];    /*    grid numbers (indexes into volt[]) on which this point 
  84.                         depends.  at least one of these will be
  85.                         negative, corresponding to a boundary value. */
  86.     int coef[4];    /*    Corresponding coefficients, scaled up by SCALE.
  87.                         The sum of these four values is approximately SCALE. */
  88.     } boundary_rule;
  89.  
  90. union     /*    these can share space */
  91.     {boundary_rule    b_rule[MAX_CROSS+3];
  92.     crossing        cross_list[MAX_CROSS+3];
  93.     } bo;
  94.  
  95. show_b(b) boundary_rule *b;
  96. {    int i, sum;
  97.     sum=0;
  98.     for (i=0; i<4; i++) sum += b->coef[i];
  99.     printf("rule %d: <>%2d>*%d + <<%2d>*%d + <^%2d>*%d + <v%2d>*%d  (ttl wt=%d%s)\n", 
  100.         b-bo.b_rule, 
  101.         b->var[0], b->coef[0], 
  102.         b->var[1], b->coef[1], 
  103.         b->var[2], b->coef[2], 
  104.         b->var[3], b->coef[3], 
  105.         sum, 
  106.         abs(sum-SCALE)>2?", WRONG!":"");
  107.     i=0;
  108. }
  109.  
  110. int lc;            /*    index of 1st unused entry in bo.cross_list    */
  111. int lb;            /*    index of 1st unused entry in bo.b_rule    */
  112.  
  113. /*
  114.     cells are numbered as follows...
  115.  
  116.                      "top"
  117.             0   1   2   3   ...   w-1                <--- ymin
  118.     "left"  w  w+1 w+2                   "right"
  119.            2w
  120.  
  121.                                ...hw-1               <--- ymax
  122.                     "bottom"
  123.             ^                      ^
  124.             |                      |
  125.            xmin                   xmax
  126.  
  127.  
  128.     cell types are as follows...
  129.  
  130.             2 3 3 3 4
  131.             9 1 1 1 5
  132.             9 1 1 1 5
  133.             9 1 1 1 5
  134.             8 7 7 7 6
  135.  
  136.     and type 0 cells have a specified potential (part of the boundary
  137.     conditions) and are never updated.  There are special cell types
  138.     are for boundaries where the normal gradient vanishes:
  139.  
  140.              10  10  10  10
  141.            13              11
  142.            13              11
  143.            13              11
  144.              12  12  12  12
  145. */
  146. int alpha=0, alphap1=1;
  147. int *type; 
  148. int *volt;            /* volt[-1]...volt[-boundaries] are boundary voltages
  149.                         volt[0]...volt[hw-1] are voltages on the grid */
  150. int debugging=0;
  151. int pass=0;
  152. int quiet=0;
  153. int show_resid=0;    /* nonzero if residuals are to be shown */
  154. int grounded=0;        /* nonzero if edge of grid is to be grounded */
  155. int show_boundaries=0;    /* nonzero if boundaries are to be shown */
  156. int show_potentials=0;    /* nonzero if calculated potentials are to be shown */
  157. int left_symmetric=0;    /* nonzero if y gradient vanishes on left edge of grid */
  158. int right_symmetric=0;    /* nonzero if y gradient vanishes on right edge of grid */
  159. int lower_symmetric=0;    /* nonzero if x gradient vanishes on lower edge of grid */
  160. int upper_symmetric=0;    /* nonzero if x gradient vanishes on upper edge of grid */
  161. int done=0;
  162. long work=0;        /* number of cell updates so far */
  163.  
  164. double xmin=1.e30;
  165. double xmax=-1.e30;
  166. double ymin=1.e30;
  167. double ymax=-1.e30;
  168. double zmin=1.e30;
  169. double zmax=-1.e30;
  170. double zavg;
  171. double margin=0.;        /* minimum amount to expand grid beyond given boundaries */
  172. double *x;                /* pointer to data area for boundaries */
  173. double *y;                /* pointer to data area for boundaries */
  174. double xscale, yscale, zscale;
  175. double edge;
  176. double aspect, correct;
  177. double p_voltage[MAXBOUNDARIES];
  178. double residual();
  179.  
  180. int boundaries=0;        /* number of boundaries in input data */
  181. int x_arguments=0;
  182. int y_arguments=0;
  183. int height=25;        /* default height of grid in cells */
  184. int width=25;        /* default width of grid in cells */
  185. int cycles=75;        /* default # relaxation cycles */
  186. int n;                /* number of entries in x and y */
  187. int index_array[MAXBOUNDARIES];    /* indexes into x and y */
  188. int *p_data=index_array;
  189. char outname[40];
  190.  
  191. FILE file;
  192. FILE ofile;
  193.  
  194.  
  195. main (argc, argv) int argc; char **argv;
  196. {    int i, j, k, hw, col;
  197.     double yval;
  198.     char *s;
  199.     read_data(argc, argv);
  200. /*
  201.     printf("there are %d points...\n", n);
  202.     k=0;
  203.     for (i=0; i<n; i++)
  204.         {printf("%f %f\n", x[i], y[i]);
  205.         if(i==p_data[k]) printf("       ...at voltage %f\n", p_voltage[k++]);
  206.         }
  207. */
  208.     if(s=index(argv[1], '.')) strncpy(outname, argv[1], s-argv[1]);
  209.     else strcpy(outname, argv[1]);
  210.     strcat(outname, ".pot");
  211.     ofile=fopen(outname, "w");
  212.     if(ofile==0) {printf("can\'t create output file\n"); exit(1);}
  213.         /* increase grid height or width to approximate aspect ratio of data */
  214.     aspect=(ymax-ymin)/(xmax-xmin);
  215.     if(height-1<(width-1)*aspect) 
  216.         height=1|(int)((width-1)*aspect+1);  /* odd number allows recursion */
  217.     else width=1|(int)((height-1)/aspect+1);
  218.         /* for periodic x boundaries, grid width must match data exactly */
  219.     if(left_symmetric && right_symmetric)
  220.         while((height-1)<(width-1)*aspect)
  221.             height += 2;
  222.         /* increase data width or height to match grid */
  223.     if(left_symmetric && right_symmetric && upper_symmetric && lower_symmetric
  224.         && fabs((height-1)-(width-1)*aspect)>.002*(height-1))
  225.         {printf("doubly periodic boundary (-xp and -yp) conditions require\n");
  226.         printf("(height-1)/(width-1) = (ymax-ymin)/(xmax-xmin)\n");
  227.         exit(1);
  228.         }
  229.     if((height-1)<(width-1)*aspect)
  230.         {correct=(ymax-ymin)*(width-1)/(height-1.)-(xmax-xmin);
  231.         if(left_symmetric) xmax += correct;
  232.         else {xmin -= correct/2.; xmax += correct/2.;}
  233.         }
  234.     else
  235.         {correct=(xmax-xmin)*(height-1)/(width-1.)-(ymax-ymin);
  236.         if(lower_symmetric) ymax += correct;
  237.         else {ymin -= correct/2.; ymax += correct/2.;}
  238.         }
  239.     hw=height*width;
  240.     if(cycles<height+width) cycles=height+width;
  241.     type=malloc(hw*sizeof(int));
  242.     volt=malloc((boundaries+hw)*sizeof(int)); 
  243.     if(!type || !volt)
  244.         {fprintf(stderr, "not enough workspace for %d by %d grid", height, width);
  245.         exit(1);
  246.         }
  247.     volt += boundaries;        /* allow space for the boundary potentials */
  248.     laplace(type, volt, width, height, cycles);
  249.     if(show_potentials) show_volt("", type, volt, width, height);
  250. #ifdef FIXED
  251.     charges(type, volt, width, height);
  252. #endif
  253.     if(!quiet) printf("writing result to file %s\n", outname);
  254.     if(-1==fprintf(ofile, "%10.4g%10.4g%5d    minimum and maximum x values, width\n", xmin, xmax, width)) exit(1);
  255.     if(-1==fprintf(ofile, "%10.4g%10.4g%5d    minimum and maximum y values, height\n", ymin, ymax, height)) exit(1);
  256.     yval=ymin;
  257.     for (i=0; i<height; i++)
  258.         {col=0;
  259.         for (j=0; j<width; j++)
  260.             {if(-1==fprintf(ofile, "%10.4g", volt[i*width+j]/zscale+zavg)) 
  261.                 exit(1);
  262.             if(++col>=7)
  263.                 {if(-1==fprintf(ofile, "\n")) exit(1);
  264.                 col=0;
  265.                 }
  266.             }
  267.         if(col)
  268.             if(-1==fprintf(ofile, "\n")) exit(1);
  269.         }
  270.     fclose(ofile);
  271. }
  272.  
  273. laplace (type, volt, width, height, cycles) int *type; int *volt, width, height, cycles;
  274. {
  275. #ifdef DESMET
  276.     long start, tics();
  277. #endif
  278.     int i, j, k, hw;
  279.     hw=height*width;
  280.     if(height<width) i=height; else i=width;
  281.     if(i>4 && (width&1)==1 && (height&1)==1)
  282.         {laplace(type, volt, width/2+1, height/2+1, cycles/2);
  283. /*    show_volt("potentials on coarse grid", type, volt, width/2+1, height/2+1); */
  284.         i=hw; j=(height/2+1)*(width/2+1);
  285.         while(i>width)            /* copy voltages into refined grid */
  286.             {i--; j--;
  287.             volt[i]=volt[i-width]=volt[j];
  288.             if(i%width==0) i -= width;
  289.             else {volt[i-1]=volt[i-width-1]=volt[j]; i--;}
  290.             }
  291.         while(i>1)
  292.             {i--; j--;
  293.             volt[i]=volt[i-1]=volt[j];
  294.             i--;
  295.             }
  296. /*    show_volt("potentials on fine grid, fully transferred", type, volt, width, height); */
  297.         }
  298.     else
  299.         {for (i=0; i<hw; i++) volt[i]=0;
  300.         }
  301.     if(!quiet) printf("pass %d: height=%d, width=%d, cycles=%d\n", 
  302.         ++pass, height, width, cycles);
  303. #ifdef DESMET
  304.     start=tics();
  305. #endif
  306.     fill(type, volt, width, height);
  307.     k=0;
  308.     if(show_boundaries)
  309. /*        {for (i=height-1; i>=0; i--) */
  310.         {for (i=0; i<height; i++)
  311.             {k=i*width;
  312.             printf("%3d: ", k);
  313.             for (j=0; j<width; j++)
  314.                 printf("%5d", type[k++]);
  315.             printf("\n");
  316.             }
  317.         if(debugging) 
  318.             {int c, g=0, t;
  319.             while(1)
  320.                 {t=type[g];
  321.                 printf("<%2d>: type %d     \n", g, t);
  322.                 if(t<0) show_b(&bo.b_rule[-t]); 
  323.                 else if(t==0) printf("potential=%f\033K\n", volt[g]/zscale+zavg);
  324.                 else printf("\033K\n");
  325.                 printf(">  \008"); c=getchar();
  326.                 if(c==56) {if(g>=width) g-=width;}                        /* up */
  327.                 else if(c==50) {if(g<height*width-width) g+=width;}        /* down */
  328.                 else if(c==54) {if(g<height*width-1) g+=1;}                /* right */
  329.                 else if(c==52) {if(g>0) g-=1;}                            /* left */
  330.                 else if(c=='q') break;                        
  331.                 printf("\015\033A\033A");
  332.                 }
  333.             }
  334.         }
  335. /*    show_volt("", type, volt, width, height); */
  336.     while(cycles--)
  337.         {if(!show_resid && !quiet)
  338.             printf("\015 %3d cycles to go  ", cycles);
  339.         relax(type, volt, width, height);
  340.         work += hw;
  341.         if(show_resid && ++done%5==0)
  342.             printf("%D   %f\n", work, residual(type, volt, width, height));
  343.         }
  344. #ifdef DESMET
  345.     if(!quiet) printf("\015    finished in %4.2f sec  \n", .01*(tics()-start));
  346. #else
  347.     if(!quiet) printf("\015    finished               \n");
  348. #endif
  349. }
  350.  
  351. relax (type, v, w, h) int *type; int *v, w, h;
  352. {    int i, wh;
  353.     int hm1, wm1, j;
  354.     boundary_rule *rule;
  355.     long sum;
  356. #ifdef DEBUG3
  357.     char *typ=type;
  358. #endif
  359. #ifdef SWITCH
  360.     i=wh=w*h;
  361.     while(i--)            /* loop over all grid points */
  362.         {if(*type>0)    /* grid point not next to boundaries */
  363.             {switch(*type)
  364.                 {
  365.                 case  1: v[0]= (  v[-w] + v[-1] +     v[1] + v[w] + 2)/4; break;
  366.                 case  2: v[0]= (                      v[1] + v[w] + 1)/2; break;
  367.                 case  3: v[0]= (        2*v[-1] + 2*v[1]   + v[w] + 2)/5; break;
  368.                 case  4: v[0]= (          v[-1]            + v[w] + 1)/2; break;
  369.                 case  5: v[0]= (2*v[-w] + v[-1]          + 2*v[w] + 2)/5; break;
  370.                 case  6: v[0]= (  v[-w] + v[-1]                   + 1)/2; break;
  371.                 case  7: v[0]= (  v[-w] + 2*v[-1] + 2*v[1]        + 2)/5; break;
  372.                 case  8: v[0]= (  v[-w]         +     v[1]        + 1)/2; break;
  373.                 case  9: v[0]= (2*v[-w]         +   v[1] + 2*v[w] + 2)/5; break;
  374.                 case 10: v[0]= (2*v[ w]   + v[-1] + v[1]           + 2)/4; break;
  375.                 case 11: v[0]= (  v[-w] + 2*v[-1]          + v[ w] + 2)/4; break;
  376.                 case 12: v[0]= (            v[-1] + v[1] + 2*v[-w] + 2)/4; break;
  377.                 case 13: v[0]= (  v[-w]         + 2*v[1]   + v[ w] + 2)/4; break;
  378.                 }
  379.             }
  380.         else if (*type<0)    /* grid point next to boundaries */
  381.             {rule=&bo.b_rule[-*type];
  382. #ifdef DEBUG3
  383. printf("\ngrid %d: ", wh-i-1); show_b(rule);
  384. show_volt("", typ, volt, w, h);
  385. #endif
  386.             sum =     volt[rule.var[0]]*(long)rule.coef[0]
  387.                     +volt[rule.var[1]]*(long)rule.coef[1]
  388.                     +volt[rule.var[2]]*(long)rule.coef[2]
  389.                     +volt[rule.var[3]]*(long)rule.coef[3];
  390.             v[0]=sum/SCALE;
  391. #ifdef DEBUG3
  392. printf("setting potential to %5.2f\n", v[0]/zscale+zavg);
  393. #endif
  394.             }
  395. #ifdef DEBUG3
  396.         if(v[0]/zscale+zavg<zmin || v[0]/zscale+zavg>zmax || (v-volt)%11==5)
  397.             {printf("\ngrid %d, of type %d\n", v-volt, type[0]);
  398.             show_star(v, w);
  399.             if(type[0]<0)  
  400.                 {int i;
  401.                 show_b(&bo.b_rule[-*type]);
  402.                 printf("\nboundary potentials:");
  403.                 for (i=1; i<=boundaries; i++)
  404.                     printf("  <%3d>=%6d=>%6.3f", 
  405.                         -i, volt[-i], volt[-i]/zscale+zavg);
  406.                 printf("\n");
  407.                 printf("<%3d> %6d * %5ld +\n", rule.var[0], volt[rule.var[0]], (long)rule.coef[0]);
  408.                 printf("<%3d> %6d * %5ld +\n", rule.var[1], volt[rule.var[1]], (long)rule.coef[1]);
  409.                 printf("<%3d> %6d * %5ld +\n", rule.var[2], volt[rule.var[2]], (long)rule.coef[2]);
  410.                 printf("<%3d> %6d * %5ld = sum = %ld => %ld => %6.3f\n", 
  411.                     rule.var[3], volt[rule.var[3]], (long)rule.coef[3],
  412.                     sum, sum/SCALE, sum/SCALE/zscale+zavg);
  413.                 }
  414.             }
  415. #endif
  416.  
  417.         v++; type++;
  418.         }
  419. #else /* not SWITCH */
  420.     hm1=h-1;
  421.     wm1=w-1;
  422.     if(*type++)
  423.         v[0]= (                      v[1] + v[w] + 1)/2;
  424.     v++;
  425.     for (j=1; j<wm1; j++)
  426.         {if(*type++)
  427.             v[0]= (           v[-1] + v[1] + v[w] + 1)/3;
  428.         v++;
  429.         }
  430.     if(*type++)
  431.         v[0]= (           v[-1]           + v[w] + 1)/2;
  432.     v++;
  433.     for (i=1; i<hm1; i++)
  434.         {
  435.         if(*type++)
  436.             v[0]= (v[-w]            + v[1] + v[w] + 1)/3;
  437.         v++;
  438.         for (j=1; j<wm1; j++)
  439.             {if(*type++)
  440.                 v[0]= (v[-w] + v[-1] + v[1] + v[w] + 2)/4;
  441.             v++;
  442.             }
  443.         if(*type++)
  444.             v[0]= (v[-w] + v[-1]           + v[w] + 1)/3;
  445.         v++;
  446.         }
  447.     if(*type++)
  448.         v[0]= (v[-w]            + v[1]           + 1)/2;
  449.     v++;
  450.     for (j=1; j<wm1; j++)
  451.         {if(*type++)
  452.             v[0]= (v[-w] + v[-1] + v[1]           + 1)/3;
  453.         v++;
  454.         }
  455.     if(*type++)
  456.         v[0]= (v[-w] + v[-1]                     + 1)/2;
  457.     v++;
  458.  
  459. #endif
  460. }
  461.  
  462. show_star(x, w) int *x, w;
  463. {    printf("                      "); show_point(x-w); printf("\n");
  464.     show_point(x-1); show_point(x); show_point(x+1); printf("\n");
  465.     printf("                      "); show_point(x+w); printf("\n");
  466. }
  467.  
  468. show_point(x) int *x;
  469. {    printf("  <%3d>=%6d=>%6.3f", x-volt, *x, *x/zscale+zavg);
  470. }
  471.  
  472. double residual (type, volt, w, h) int *type; int *volt, w, h;
  473. {    int i, hw, n;
  474.     int hm1, wm1, j;
  475.     long r, d;
  476.  
  477.     n=0;
  478.     i=hw=w*h;
  479.     r=0L;
  480.     while(i--)
  481.         {switch(*type++)
  482.             {case 0: n++; break;
  483.             case 1: d=volt[0]-(volt[-w] + volt[-1] + volt[1] + volt[w] + 2)/4; r+=d*d; break;
  484.             case 2: d=volt[0]-(                      volt[1] + volt[w] + 1)/2; r+=d*d; break;
  485.             case 3: d=volt[0]-(           volt[-1] + volt[1] + volt[w] + 1)/3; r+=d*d; break;
  486.             case 4: d=volt[0]-(           volt[-1]           + volt[w] + 1)/2; r+=d*d; break;
  487.             case 5: d=volt[0]-(volt[-w] + volt[-1]           + volt[w] + 1)/3; r+=d*d; break;
  488.             case 6: d=volt[0]-(volt[-w] + volt[-1]                     + 1)/2; r+=d*d; break;
  489.             case 7: d=volt[0]-(volt[-w] + volt[-1] + volt[1]           + 1)/3; r+=d*d; break;
  490.             case 8: d=volt[0]-(volt[-w]            + volt[1]           + 1)/2; r+=d*d; break;
  491.             case 9: d=volt[0]-(volt[-w]            + volt[1] + volt[w] + 1)/3; r+=d*d; break;
  492.             }
  493.         volt++;
  494.         }
  495.     if (hw<=n)
  496.         return 1.;
  497.     return sqrt(((double)r)/(hw-n));
  498. }
  499.  
  500. #ifdef FIXED
  501. /*
  502.     this no longer works...could be changed to solve linear system for
  503.     the gradient at center of each grid point next to a boundary and
  504.     integrate, but it's easier to let CONTOUR calculate the integral.
  505. */
  506. charges (type, volt, width, height) int *type; int *volt, width, height;
  507. {    int v, i, plate=0, out=0, hw;
  508.     long charge;
  509.     double C, Z, q[2], plate_voltage[2], sum=0.;
  510.     hw=height*width;
  511.     printf("\nComputing surface integral of grad V over each boundary\n");
  512.     printf("(integral of -grad V dot N, where N is normal to the boundary)...\n");
  513.     printf("\nboundary    potential    surface integral     charge\n");
  514.     while(1)
  515.         {for (i=0; i<hw; i++)
  516.             if(type[i]==0)
  517.                 break;
  518.         if(i==hw) break;
  519.         v=volt[i];
  520.         plate_voltage[plate&1]=v/zscale+zavg;
  521.         charge=0.;
  522.         if(i%width) charge += v-volt[i-1];
  523.         if(i>=width) charge += v-volt[i-width];
  524.         if(i<hw-width) charge += v-volt[i+width];
  525.         if(i%width != width-1) charge += v-volt[i+1];
  526.         type[i]=23;
  527.         for (i++; i<hw; i++)
  528.             {if(type[i]==0 && v==volt[i])
  529.                 {if(i%width) charge += v-volt[i-1];
  530.                 if(i>=width) charge += v-volt[i-width];
  531.                 if(i<hw-width) charge += v-volt[i+width];
  532.                 if(i%width != width-1) charge += v-volt[i+1];
  533.                 type[i]=23;        /* don't look at this cell again */
  534.                 }
  535.             }
  536.         sum += (q[plate&1]=(charge/zscale)*PERMITTIVITY);
  537.         printf("%5d %14.2g %15.3g %15.3g\n", 
  538.             plate, v/zscale+zavg, q[plate&1]/PERMITTIVITY, q[plate&1]);
  539.         plate++;
  540.         }
  541.     printf("    sum (should be zero)%12.3g %15.3g\n\n", sum/PERMITTIVITY, sum);
  542.     if(plate==2)
  543.         {q[0]=fabs(q[0]);
  544.         q[1]=fabs(q[1]);
  545.         C=MAX(q[1], q[0])/fabs(plate_voltage[1]-plate_voltage[0]);
  546.         Z=1./C/2.9979e8;
  547.         printf("characteristic impedance is %5.2f ohms\n", Z);
  548.         printf("capacitance is %7.3g farad/m\n", C);
  549.         printf("inductance is %7.3g henries/m\n", Z/2.9979e8);
  550.         }
  551. }
  552. #endif
  553. /*
  554.     the equations relating a grid point and its four neighbors, in the
  555.     absence of a boundary, are:
  556.     
  557.         ( 0  1  1  1) (Uy   )   (Uright)
  558.         ( 0 -1  1  1) (Ux   ) = (Uleft )
  559.          ( 1  0 -1  1) (Uxx  )   (Uup   )
  560.         (-1  0 -1  1) (Uhere)   (Udown )
  561.  
  562.     where    Ux is h*(d/dx U)
  563.             Uy is h*(d/dy U)
  564.             Uxx is .5*h**2*(d^2/dx^2 U))
  565.  
  566.     near a boundary, the equations are:
  567.  
  568.         ( 0  R  R*R  1) (Uy   )   (Uright)
  569.         ( 0 -L  L*L  1) (Ux   ) = (Uleft )
  570.          ( U  0 -U*U  1) (Uxx  )   (Uup   )
  571.         (-D  0 -D*D  1) (Uhere)   (Udown )
  572.  
  573.     where neighbors to right, left, up, and down are at distances R, L, 
  574.     U, and D respectively.
  575.  
  576.     FILL inverts the matrix and saves the last row of the inverse.  The
  577.     dot product of that row with the right hand side of the above
  578.     equation yields Uhere, the revised value of the potential at the
  579.     center.
  580.     
  581. */
  582.  
  583. double lhs[16]=
  584.     {0., 1., 1., 1., 
  585.      0., -1., 1., 1., 
  586.      1., 0., -1., 1., 
  587.     -1., 0., -1., 1.};
  588.  
  589. #define LEFT_SIDE(g)    ((g)%w==0)
  590. #define RIGHT_SIDE(g)    ((g)%w==w-1)
  591. #define TOP_SIDE(g)        ((g)<w)
  592. #define LOWER_SIDE(g)    ((g)>=hw-w)
  593.  
  594.     /*    fill type and voltage arrays based on user's boundary conditions */
  595. fill (type, volt, w, h) int *type; int *volt, w, h;
  596. {    boundary_rule *ru;
  597.     crossing *this, *last;
  598.     double invert(), decomp();
  599.     double a[4][4], ai[4][4], cond, condp1, *p, f;
  600.     double xt, yt, xt2, yt2;
  601.     int i, j, g, zero, *ip;
  602.     long weight; 
  603.     int compare();
  604.     int i, k, hw, zz;
  605.  
  606.     if(!quiet) printf("    following boundaries \n");
  607.     hw=h*w;
  608.     xscale=(w-1)/(xmax-xmin);
  609.     yscale=(h-1)/(ymax-ymin);    /* same as xscale */
  610.     if(fabs(xscale-yscale)>.001*fabs(yscale))
  611.         {printf("grid isn't square... xmin=%f, xmax=%f, w=%d, xscale=%f\n", 
  612.                 xmin, xmax, w, xscale);
  613.         printf("                     ymin=%f, ymax=%f, w=%d, yscale=%f\n", 
  614.                 ymin, ymax, h, yscale);
  615.         }
  616.     edge=1./xscale;
  617. #ifdef DEBUG1
  618.     printf("xscale=%f, yscale=%f, edge=%f\n", xscale, yscale, edge); 
  619. #else
  620. #ifdef DEBUG2
  621.     printf("xscale=%f, yscale=%f, edge=%f\n", xscale, yscale, edge); 
  622. #endif /* DEBUG2 */
  623. #endif /* DEBUG1 */
  624.  
  625.     for (i=0; i<hw; i++)    type[i]=1;        /* interior */
  626.     if(grounded)
  627.         {zero=floor((0.-zavg)*zscale+.499);
  628.         for (i=w; i<hw; i+=w)    {type[i]=0; volt[i]=zero;}    /* left */
  629.         for (i=0; i<w; i++)        {type[i]=0; volt[i]=zero;}    /* top */
  630.         for (i=w-1; i<hw; i+=w)    {type[i]=0; volt[i]=zero;}    /* right */
  631.         for (i=hw-w; i<hw; i++)    {type[i]=0; volt[i]=zero;}    /* bottom */
  632.         }
  633.     else
  634.         {if(left_symmetric)
  635.             for (i=w; i<hw; i+=w)    type[i]=13;        /* left */
  636.         else
  637.             for (i=w; i<hw; i+=w)    type[i]=9;        /* left */
  638.         for (i=0; i<w; i++)        type[i]=3;            /* top */
  639.         for (i=w-1; i<hw; i+=w)    type[i]=5;            /* right */
  640.         if(lower_symmetric)
  641.             for (i=hw-w; i<hw; i++)    type[i]=12;        /* bottom */
  642.         else
  643.             for (i=hw-w; i<hw; i++)    type[i]=7;        /* bottom */
  644.         type[0]=2; type[w-1]=4; type[hw-w]=8; type[hw-1]=6;        /* corners */
  645.         }
  646.     i=k=0;
  647.     lc=2;    /* bo.cross_list is empty */
  648.     while(i<n)
  649.         {xt=x[i]-xmin;
  650.         yt=y[i]-ymin;
  651.         zz=floor((p_voltage[k]-zavg)*zscale+.499);
  652.         volt[-k-1]=zz;
  653.         do    {xt2=x[i]-xmin;
  654.             yt2=y[i]-ymin;
  655.             mark(w, h, xt, yt, xt2, yt2, -k-1);
  656. /*
  657.             printf("%f, %f -> ", (x[i]-xmin)*xscale, (y[i]-ymin)*yscale);
  658.             printf("mark(%d, %d, %d, %d, %d, %d, %d)\n", w, h, xx, yy, xx2, yy2, zz); 
  659. */
  660.             xt=xt2; yt=yt2; 
  661.             i++;
  662.             } while(i<=p_data[k]);
  663.         k++;
  664.         }
  665.     if(!quiet) printf("    boundary crosses grid %d times\n", lc-1);
  666. /*----------------  sort list  -------------*/
  667.     hsort(&bo.cross_list[2], lc-2, sizeof(crossing), &compare);
  668. /*----------------  solve each problem  -------------*/
  669. #ifdef DEBUG2
  670.     printf("crossings after sorting\n");
  671.     for (i=2; i<17; i++) show_c(&bo.cross_list[i], w); getchar();
  672. #endif /* DEBUG2 */
  673. #ifdef DEBUG4
  674.         {int g; double xx, yy, d; FILE *out;
  675.         out=fopen("points","w");
  676.         for (i=2; i<=lc; i++) 
  677.             {g=bo.cross_list[i].grid;
  678.             xx=(g%w)*edge+xmin; yy=(g/w)*edge+ymin;
  679.             d=bo.cross_list[i].dist*edge;
  680.             switch(bo.cross_list[i].direct)
  681.                 {case 0: xx+=d;    break;    /* right */
  682.                 case 1: xx-=d;    break;    /* left */
  683.                 case 2: yy-=d;    break;    /* above */
  684.                 case 3: yy+=d;    break;    /* below */
  685.                 }
  686.             fprintf(out, "%10.3f %10.3f\n", xx, yy);
  687.             }
  688.         fclose(out);
  689.         }
  690. #endif /* DEBUG4 */
  691.  
  692.     lb=1;    /* bo.b_rule is empty */
  693.     this=&bo.cross_list[2];
  694.     last=&bo.cross_list[lc];
  695.     last->grid=hw;            /* sentinal */
  696.     while(this<last)
  697.         {p=a;
  698.         for (i=0; i<16; i++) p[i]=lhs[i];
  699.         g=this->grid;
  700.         if(g<0) {this++; continue;}
  701.         if(type[g]==0) {while((++this)->grid==g){} continue;}
  702.         if(g>=hw) break;
  703. #ifdef DEBUG2
  704. printf("rule %d, boundary conditions for grid point %d\n", lb, g);
  705. #endif
  706.                 /* set voltage for this grid to that of the 
  707.                     first boundary in the list (the closest one)    */
  708.         volt[g]=floor((p_voltage[-this->volt-1]-zavg)*zscale+.499);
  709.         ru=&bo.b_rule[lb];
  710.                     /* initialize variables to refer to neighbors */
  711.         ru->var[0]=g+1;        /* right */
  712.         ru->var[1]=g-1;        /* left */
  713.         ru->var[2]=g-w;        /* up */
  714.         ru->var[3]=g+w;        /* down */
  715.         do    {j=this->direct;
  716. #ifdef DEBUG2
  717. printf("    applying:"); show_c(this, w);
  718. #endif
  719.             f=this->dist;
  720.             if(f<edge*.01)      /* if boundary is very close, */
  721.                 {type[g]=0;     /* declare the point to be type 0 */
  722.                 while((++this)->grid==g){} 
  723.                 goto NEXT_GRID;
  724.                 }
  725.             ru->var[j]=this->volt;    /* point to boundary instead of neighbor */
  726.             a[j][0] *= f;
  727.             a[j][1] *= f;
  728.             a[j][2] *= f*f;
  729.             while((++this)->grid==g && this->direct==j) 
  730.                 {
  731. #ifdef DEBUG2
  732. printf("    skipping:"); show_c(this, w);
  733. #endif
  734.                 }
  735.             } while (this->grid==g);
  736. /*        show_m("boundary condition matrix", a); getchar();    */
  737.         cond=invert(4, 4, a, ai);
  738. /*        printf("condition number is %f\n", cond);
  739.         show_m("matrix inverse", ai);    */
  740.         condp1=cond+1.;
  741.         if(cond==condp1)
  742.             {fprintf(stderr, 
  743.                 "can\'t solve linear system for boundary conditions\n");
  744.             while((--this)->grid==g) {fprintf(stderr, "        "); show_c(this, w);}
  745.             exit(1);
  746.             }
  747.         ru->coef[0]=ai[3][0]*SCALE+.5;  /* unrolled for speed */
  748.         ru->coef[1]=ai[3][1]*SCALE+.5;
  749.         ru->coef[2]=ai[3][2]*SCALE+.5;
  750.         ru->coef[3]=ai[3][3]*SCALE+.5;
  751.         if(LOWER_SIDE(g)||LEFT_SIDE(g)||TOP_SIDE(g)||RIGHT_SIDE(g))
  752.             {
  753. #ifdef DEBUG2
  754.             printf("rule for grid point %d, ", g);
  755.             if(LOWER_SIDE(g)) printf("LOWER ");
  756.             if(LEFT_SIDE(g)) printf("LEFT ");
  757.             if(TOP_SIDE(g)) printf("TOP ");
  758.             if(RIGHT_SIDE(g)) printf("RIGHT ");
  759.             printf("\n");
  760.             printf("before adjustment ");show_b(ru);
  761. #endif
  762.             }
  763.         if(LOWER_SIDE(g))
  764.             {if(lower_symmetric)
  765.                 ru->coef[2] += ru->coef[3];
  766.             else
  767.                 {weight=ru->coef[2]=(ru->coef[2] + ru->coef[3])/4;
  768.                 weight = (SCALE*(long)SCALE)/(weight+ru->coef[0]+ru->coef[1]);
  769.                 ru->coef[0] = (ru->coef[0]*weight)/SCALE;
  770.                 ru->coef[1] = (ru->coef[1]*weight)/SCALE;
  771.                 ru->coef[2] = (ru->coef[2]*weight)/SCALE;
  772.                 }
  773.             ru->coef[3]=0;
  774.             }
  775.         else if(TOP_SIDE(g))
  776.             {if(upper_symmetric)
  777.                 ru->coef[3] += ru->coef[2];
  778.             else
  779.                 {weight=ru->coef[3]=(ru->coef[2] + ru->coef[3])/4;
  780.                 weight = (SCALE*(long)SCALE)/(weight+ru->coef[0]+ru->coef[1]);
  781.                 ru->coef[0] = (ru->coef[0]*weight)/SCALE;
  782.                 ru->coef[1] = (ru->coef[1]*weight)/SCALE;
  783.                 ru->coef[3] = (ru->coef[3]*weight)/SCALE;
  784.                 }
  785.             ru->coef[2]=0;
  786.             }
  787.         if(LEFT_SIDE(g))
  788.             {if(left_symmetric)
  789.                 ru->coef[0] += ru->coef[1];
  790.             else
  791.                 {weight=ru->coef[0]=(ru->coef[0] + ru->coef[1])/4;
  792.                 weight = (SCALE*(long)SCALE)/(weight+ru->coef[2]+ru->coef[3]);
  793.                 ru->coef[1] = (ru->coef[1]*weight)/SCALE;
  794.                 ru->coef[2] = (ru->coef[2]*weight)/SCALE;
  795.                 ru->coef[3] = (ru->coef[3]*weight)/SCALE;
  796.                 }
  797.             ru->coef[1]=0;
  798.             }
  799.         else if(RIGHT_SIDE(g))
  800.             {if(right_symmetric)
  801.                 ru->coef[1] += ru->coef[0];
  802.             else
  803.                 {weight=ru->coef[1]=(ru->coef[1] + ru->coef[0])/4;
  804.                 weight = (SCALE*(long)SCALE)/(weight+ru->coef[2]+ru->coef[3]);
  805.                 ru->coef[1] = (ru->coef[1]*weight)/SCALE;
  806.                 ru->coef[2] = (ru->coef[2]*weight)/SCALE;
  807.                 ru->coef[3] = (ru->coef[3]*weight)/SCALE;
  808.                 }
  809.             ru->coef[0]=0;
  810.             }
  811. #ifdef DEBUG2
  812. printf("    ");show_b(ru);
  813. #endif
  814.         type[g]=-lb;
  815.         lb++;
  816. NEXT_GRID:
  817.         ;
  818.         }
  819.     if(!quiet) printf("    %d grid points have special equations\n", lb-1);
  820. }
  821.  
  822. compare(a, b) crossing *a, *b;
  823. {    int n;
  824.     if(n=(a->grid - b->grid)) return n;
  825.     if(n=(a->direct - b->direct)) return n;
  826.     if(a->dist < b->dist) return -1;
  827.     return 1;    /* if same distance away, declare a > b */
  828. }
  829.  
  830. mark (width, height, x1, y1, x2, y2, v) 
  831. int width, height,         /* dimensions of grid */
  832. v;                        /* potential of boundary segment is volt[-v-1] */
  833. double x1, y1, x2, y2;    /* beginning and end points of boundary segment */
  834. {    int q, dx, dy, ax, ay, decide, up, over, i;
  835.     int nx, nx2, ny, ny2;
  836.     double del, xe, ye, f, t;
  837. /*
  838.     if(x1<0||x1>=width||x2<0||x2>=width||y1<0||y1>=height||y2<0||y2>=height)
  839.         {fprintf(stderr, "calling sequence error: mark(%d, %d, %d, %d, %d, %d)\n", 
  840.             width, x1, y1, x2, y2, v);
  841.         fprintf(stderr, "(point outside 0<=x<%d or 0<=y<%d)\n", width, height);
  842.         exit(1);
  843.         }
  844. */
  845. #ifdef DEBUG1
  846. printf("line from  %5.2f, %5.2f  to  %5.2f, %5.2f at potential %5.2f\n", 
  847. x1, y1, x2, y2, p_voltage[-v-1]);
  848. printf("     slope = dy/dx = %5.2f/%5.2f\n", y2-y1, x2-x1);
  849. #endif
  850.             /* test whether endpoint is near a gridpoint */
  851.     nx2=x2*xscale+.5;
  852.     ny2=y2*yscale+.5;
  853.     if(fabs(nx2*edge-x2)+fabs(ny2*edge-y2)<.01*edge)
  854.         register_crossing(width, nx2, ny2, 0., 2, v);
  855.     if(x2<x1) {t=x1; x1=x2; x2=t; t=y1; y1=y2; y2=t;}
  856.             /* now, x1 <= x2 */
  857.     nx2=x2*xscale;
  858.     nx=x1*xscale+1;
  859.     if(nx<=nx2)
  860.         {xe=nx*edge;
  861.         del=yscale/(x2-x1);
  862.         ye=(y1*(x2-xe) + y2*(xe-x1))*del;
  863.         del *= (y2-y1)*edge;
  864.         for ( ; nx<=nx2; nx++)
  865.             {ny=ye;
  866.             f=ye-ny;
  867. #ifdef DEBUG1
  868. printf("    crosses y grid line at  %5.3f, %5.3f: %5.3f below <%d> = (%d,%d)\n", 
  869.     xe, ye, f, nx+ny*width, nx, ny);
  870. #endif
  871.             register_crossing(width, nx, ny, f, 3, v);        /* below this point */
  872.             register_crossing(width, nx, ny+1, 1.-f, 2, v);    /* above this point */
  873.             xe += edge;
  874.             ye += del;
  875.             }
  876.         }
  877.     if(y2<y1) {t=x1; x1=x2; x2=t; t=y1; y1=y2; y2=t;}
  878.             /* now, y1 <= y2 */
  879.     ny2=y2*yscale;
  880.     ny=y1*yscale+1;
  881.     if(ny<=ny2)
  882.         {ye=ny*edge;
  883.         del=xscale/(y2-y1);
  884.         xe=(x1*(y2-ye) + x2*(ye-y1))*del;
  885.         del *= (x2-x1)*edge;
  886.         for (  ; ny<=ny2; ny++)
  887.             {nx=xe;
  888.             f=xe-nx;
  889. #ifdef DEBUG1
  890. printf("    crosses x grid line at  %5.3f, %5.3f: %5.3f right of <%d> = (%d,%d)\n", 
  891.     xe, ye, f, nx+ny*width, nx, ny);
  892. #endif
  893.             register_crossing(width, nx, ny, f, 0, v);        /* right of this */
  894.             register_crossing(width, nx+1, ny, 1.-f, 1, v); /* left of this */
  895.             ye += edge;
  896.             xe += del;
  897.             }
  898.         }
  899. }
  900.  
  901. register_crossing(width, nx, ny, f, direction, v)
  902. int width, nx, ny, direction, v;
  903. double f;
  904. {    int g;
  905.     crossing *new;
  906.  
  907.     g=nx+ny*width;
  908.     if(type[g]==0) return;
  909.     if(f<.01)
  910.         {type[g]=0;
  911.         volt[g]=floor((p_voltage[-v-1]-zavg)*zscale+.499);
  912. #ifdef DEBUG1
  913.         printf("        <%d> set to %3.2f, type 0\n", g, p_voltage[-v-1]);
  914. #endif
  915.         return;
  916.         }
  917.     if(f>.99) return;
  918.     if(lc>=MAX_CROSS) 
  919.         {fprintf(stderr, "too many boundary crosspoints\n"); 
  920.         exit(1);
  921.         }
  922.     new=&bo.cross_list[lc++];
  923.     new->grid=g;
  924.     new->direct=direction;
  925.     new->dist=f;
  926.     new->volt=v;
  927. #ifdef DEBUG1
  928.     printf("        "); show_c(new, width);
  929. #endif
  930. }
  931.  
  932. show_volt (msg, type, volt, w, h) char *msg; int *type; int *volt, w, h;
  933. {    int i, j, k;
  934.     k=0;
  935. printf("%s\nboundaries...", msg);
  936. for (i=1; i<=boundaries; i++) printf("%d: %5.2f   ", i, volt[-i]/zscale+zavg);
  937. printf("\n");
  938. /*    printf("(press any key to abort printout)\n");
  939.     for (i=0; i<h && !csts(); i++)
  940. */
  941.     for (i=0; i<h; i++)
  942.         {k=i*w;
  943.         printf("%3d: ", k);
  944.         for (j=0; j<w; j++)
  945.             {printf("%6.2f", volt[k++]/zscale+zavg);
  946.             /* printf("%6d", volt[k++]); */
  947.             }
  948.         printf("\n");
  949.         }
  950.     for (i=0; i<h; i++)
  951.         {k=i*w;
  952.         printf("%3d: ", k);
  953.         for (j=0; j<w; j++)
  954.             printf("%6d", type[k++]);
  955.         printf("\n");
  956.         }
  957.     if(debugging) getchar();
  958. }
  959.  
  960. read_data(argc, argv) int argc; char **argv;
  961. {    int i, j, length;
  962.     double xx, yy, d, *pd, sum;
  963.     char *s, *t;
  964. #define BUFSIZE 200
  965.     static char buf[BUFSIZE];
  966.  
  967.     x=malloc(ENTRIES*sizeof(double));
  968.     y=malloc(ENTRIES*sizeof(double));
  969.     if(x==0 || y==0) {fprintf(stderr, "can\'t allocate buffer"); exit();}
  970.     if(argc>1 && streq(argv[1], "?")) help();
  971.     if(argc<=1 || *argv[1]=='-') file=stdin;
  972.     else
  973.         {if(argc>1)
  974.             {file=fopen(argv[1], "r");
  975.             if(file==0) {printf("file %s not found\n", argv[1]); exit();}
  976.             argc--; argv++;
  977.             }
  978.         else help();
  979.         }
  980.     argc--; argv++;
  981.     while(argc>0)
  982.         {i=get_parameter(argc, argv);
  983.         argv=argv+i; argc=argc-i;
  984.         }
  985.     p_data[0]=-1;
  986.     if(height<2||width<2||cycles<2)
  987.         {fprintf(stderr, "parameter too small: height=%d, width=%d, cycles=%d", 
  988.             height, width, cycles);
  989.         exit(1);
  990.         }
  991.     i=0;
  992.     while(i<ENTRIES)
  993.         {if(fgets(buf, BUFSIZE, file)==0) break;
  994.         t=buf;
  995.         while(*t && isspace(*t)) t++;
  996.         if(*t == 0) continue;        /* skip blank lines */
  997.         buf[strlen(buf)-1]=0;        /* zero out the line feed */
  998.         if(buf[0]==';') {printf("%s\n", buf); continue;}  /* skip comment */
  999.         sscanf(buf, "%F %F", &x[i], &y[i]);
  1000.         if (x[i]<xmin) xmin=x[i];
  1001.         if (x[i]>xmax) xmax=x[i];
  1002.         if (y[i]<ymin) ymin=y[i];
  1003.         if (y[i]>ymax) ymax=y[i];
  1004.         s=buf;                      /* start looking for label */
  1005.         while(*s==' ')s++;            /* skip first number */
  1006.         while(*s && (*s!=' '))s++;
  1007.         while(*s==' ')s++;            /* skip second number */
  1008.         while(*s && (*s!=' '))s++;
  1009.         while(*s==' ')s++;
  1010.         if((length=strlen(s))&&(boundaries<MAXBOUNDARIES))
  1011.             {if(*s=='\"') s++;           /* label in quotes */
  1012.             sscanf(s, "%F", &p_voltage[boundaries]);
  1013.             if(p_voltage[boundaries]<zmin) zmin=p_voltage[boundaries];
  1014.             if(p_voltage[boundaries]>zmax) zmax=p_voltage[boundaries];
  1015.             p_data[boundaries++]=i;
  1016.             }
  1017.         i++;
  1018.         }
  1019.     n=i;
  1020.     if(grounded)
  1021.         {if(0.<zmin) zmin=0.;
  1022.         if(0.>zmax) zmax=0.;
  1023.         }
  1024.     if(!boundaries || p_data[boundaries-1]!=n-1)
  1025.         {p_data[boundaries]=i-1;
  1026.         p_voltage[boundaries++]=0.;
  1027.         }
  1028.     zscale=65536./5./(zmax-zmin);
  1029.     zavg=(zmax+zmin)/2.;
  1030.     xmin -= margin;
  1031.     xmax += margin;
  1032.     ymin -= margin;
  1033.     ymax += margin;
  1034. }
  1035.  
  1036.  
  1037. /* get_parameter - process one command line option
  1038.         (return # parameters used) */
  1039. get_parameter(argc, argv) int argc; char **argv;
  1040. {    int i;
  1041.     double temp;
  1042.  
  1043.     if(streq(*argv, "-d")) {debugging=1; return 1;}
  1044.     else if(streq(*argv, "-q")) {quiet=1; return 1;}
  1045.     else if(streq(*argv, "-h"))
  1046.         {if((argc>1) && numeric(argv[1])) height=atoi(argv[1]);
  1047.         return 2;
  1048.         }
  1049.     else if(streq(*argv, "-w"))
  1050.         {if((argc>1) && numeric(argv[1])) width=atoi(argv[1]);
  1051.         return 2;
  1052.         }
  1053.     else if(streq(*argv, "-g"))
  1054.         {grounded++;
  1055.         return 1;
  1056.         }
  1057.     else if(streq(*argv, "-m"))
  1058.         {i=get_double(argc, argv, 1, &margin, &margin, &margin);
  1059.         if(margin<0.) margin=0.;
  1060.         return i;
  1061.         }
  1062.     else if(streq(*argv, "-n"))
  1063.         {if((argc>1) && numeric(argv[1])) cycles=atoi(argv[1]);
  1064.         return 2;
  1065.         }
  1066.     else if(streq(*argv, "-r")) {show_resid=1; return 1;}
  1067.     else if(streq(*argv, "-b")) {show_boundaries=1; return 1;}
  1068.     else if(streq(*argv, "-p")) {show_potentials=1; return 1;}
  1069. /*    else if(streq(*argv, "-a"))
  1070.         {i=get_double(argc, argv, 1, &temp, &temp, &temp);
  1071.         alpha=temp;
  1072.         alphap1=alpha+1;
  1073.         return i;
  1074.         }
  1075. */
  1076.     else if(streq(*argv, "-x"))
  1077.         {i=get_double(argc, argv, 2, &xmin, &xmax, &xmax);
  1078.         x_arguments=i-1;
  1079.         return i;
  1080.         }
  1081.     else if(streq(*argv, "-y"))
  1082.         {i=get_double(argc, argv, 2, &ymin, &ymax, &ymax);
  1083.         y_arguments=i-1;
  1084.         return i;
  1085.         }
  1086.     else if(streq(*argv, "-ys")) {lower_symmetric=1; return 1;}
  1087.     else if(streq(*argv, "-xs")) {left_symmetric=1; return 1;}
  1088.     else if(streq(*argv, "-xp")) {right_symmetric=left_symmetric=1; return 1;}
  1089.     else if(streq(*argv, "-yp")) {upper_symmetric=lower_symmetric=1; return 1;}
  1090.     else gripe(argv);
  1091. }
  1092.  
  1093. get_double(argc, argv, permitted, a, b, c)
  1094. int argc, permitted; char **argv; double *a, *b, *c;
  1095. {    int i=1;
  1096.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *a=atof(argv[i++]);
  1097.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *b=atof(argv[i++]);
  1098.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *c=atof(argv[i++]);
  1099.     return i;
  1100. }
  1101.  
  1102. gripe_arg(s) char *s;
  1103. {    fprintf(stderr, "argument missing for switch %s", s);
  1104.     help();
  1105. }
  1106.  
  1107. gripe(argv) char **argv;
  1108. {    fprintf(stderr, *argv); fprintf(stderr, " isn\'t a legal argument \n\n");
  1109.     help();
  1110. }
  1111.  
  1112. help()
  1113. {    printf("laplace   version %s", VERSION);
  1114.     printf("\nusage: laplace [file][options]\n");
  1115.     printf("options are:\n");
  1116. /*    printf("  -a  num            acceleration factor (default %d) \n", alpha); */
  1117.     printf("  -b               display boundaries\n");
  1118.     printf("  -g               edge of grid is grounded (potential 0)\n");
  1119.     printf("  -p               display calculated potentials\n");
  1120.     printf("  -q               quiet: no comments during calculation\n");
  1121.     printf("  -r               calculate and display residuals\n");
  1122.     printf("  -m  margin       distance to expand grid beyond given boundaries (default 0) \n");
  1123.     printf("  -n  num          number of relaxation cycles (default %d) \n", cycles);
  1124.     printf("  -h  num          height of grid in cells (default %d)\n", height);
  1125.     printf("  -w  num          width of grid in cells (default %d)\n", width);
  1126.     printf("  -x  [min [max]]  specify x range \n");
  1127.     printf("  -y  [min [max]]  specify y range \n");
  1128.     printf("  -xs              symmetric across the line x=xmin \n");
  1129.     printf("  -ys              symmetric across the line y=ymin \n");
  1130.     printf("  -xp    periodic: symmetric across lines x=xmin and x=xmax \n");
  1131.     printf("  -yp    periodic: symmetric across lines y=ymin and y=ymax \n");
  1132.     exit();
  1133. }
  1134.  
  1135. numeric(s) char *s;
  1136. {    char c;
  1137.     while(c=*s++)
  1138.         {if((c<='9' && c>='0') || c=='+' || c=='-' || c=='.') continue;
  1139.         return 0;
  1140.         }
  1141.     return 1;
  1142. }
  1143.  
  1144. #ifdef DESMET
  1145.  
  1146. /*
  1147.     tics - report hundreths of a second since midnight
  1148. */
  1149.  
  1150. long tics()
  1151. {    int hr, min, sec, hun;
  1152.  
  1153.     _timer(&hr, &min, &sec, &hun);
  1154.     return (( (long) hr*60+min)*60+sec)*100+hun;
  1155. }
  1156.  
  1157. _timer()
  1158. {
  1159. #asm
  1160.     mov    ah, 2ch
  1161.     int    21h
  1162.     mov    bx, [bp+4]
  1163.     mov    [bx], ch        ;hours
  1164.     mov    byte [bx+1], 0
  1165.     mov    bx, [bp+6]
  1166.     mov    [bx], cl        ;minutes
  1167.     mov    byte [bx+1], 0
  1168.     mov    bx, [bp+8]
  1169.     mov    [bx], dh        ;seconds
  1170.     mov    byte [bx+1], 0
  1171.     mov    bx, [bp+10]
  1172.     mov    [bx], dl        ;hundreths
  1173.     mov    byte [bx+1], 0
  1174. #endasm
  1175. }
  1176.  
  1177. #endif
  1178.  
  1179. show_m(s, a) char *s; double *a;
  1180. {    int i, j;
  1181.     printf("%s\n", s);
  1182.     for (i=0; i<4; i++)
  1183.         {for (j=0; j<4; j++) printf("%8.2f ", *a++);
  1184.         printf("\n");
  1185.         }
  1186. }
  1187.  
  1188.  
  1189. char *c_label[4]={"right of", "left of", "above", "below"};
  1190. show_c(c, w) crossing *c; int w;
  1191. {    double xx, yy, d;
  1192.     int g;
  1193.     g=c->grid;
  1194.     xx=(g%w)*edge+xmin; yy=(g/w)*edge+ymin;
  1195.     d=c->dist*edge;
  1196.     switch(c->direct)
  1197.         {case 0: xx+=d;    break;    /* right */
  1198.         case 1: xx-=d;    break;    /* left */
  1199.         case 2: yy-=d;    break;    /* above */
  1200.         case 3: yy+=d;    break;    /* below */
  1201.         }
  1202.     printf("boundary value %d is %4.3f %s grid %d =(%4.3f,%4.3f)",
  1203.     c->volt, c->dist, c_label[c->direct], c->grid, xx, yy);
  1204. }
  1205.  
  1206.  
  1207.